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Abstract 

We present new splitting methods designed for the numerical integra- 
tion of ncar-intcgrablc Hamiltonian systems, and in particular for plan- 
etary N-body problems, when one is interested in very accurate results 
over a large time span. We derive in a systematic way an independent 
set of necessary and sufficient conditions to be satisfied by the coefficients 
of splitting methods to achieve a prescribed order of accuracy. Splitting 
methods satisfying such (generalized) order conditions are appropriate in 
particular for the numerical simulation of the Solar System described in 
Jacobi coordinates. We show that, when using Poincare Heliocentric co- 
ordinates, the same order of accuracy may be obtained by imposing an 
additional polynomial equation on the coefficients of the splitting method. 
We construct several splitting methods for each of the two sets of coordi- 
nates by solving the corresponding systems of polynomial equations and 
finding the optimal solutions. The experiments reported here indicate that 
the efficiency of our new schemes in both sets of canonical coordinates is 
similar, and clearly superior to previous integrators when high accuracy is 
required. 
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1 Introduction 



Symplectic integrators have several features that turn out to be particularly 
appropriate when integrating numerically for long times evolution problems in 
dynamical astronomy. They preserve by construction the symplectic structure 
of the original Hamiltonian problem, so that the numerical solution inherits the 
qualitative properties of the exact one |18| . In particular, by using backward 
error analysis, it is possible to prove that this numerical solution is in fact 
exponentially close to the exact solution of a modified Hamiltonian. Moreover, 
although the energy is not conserved along the trajectory, the error introduced 
by a symplectic method of order r used with constant step size r is of order 
0(r r ) for exponentially long time intervals under rather general assumptions, 
whereas the error in position typically grows linearly with time [6]. 

Assume that, as is often the case, the Hamiltonian function is of the form 
H(q,p) = T(p) + U(q), where the potential energy U(q) depends on positions 
and the kinetic energy T(p) is a function of the conjugate momenta. Then 
the equations of motion corresponding to T{p) are trivially solvable, and the 
same happens with U(q). By composing the flows of these two special Hamil- 
tonian systems one gets a symplectic first order approximation to the exact 
flow. This simple composition constitutes an example of a symplectic splitting 
method. Higher order approximations can be obtained by composing the flows 
corresponding to T(p) and U(q) with certain coefficients obtained by solving 
the so-called order conditions [H]. There exist in the literature a vast num- 
ber of high order integrators constructed along this line (see, e.g., [I], [6], and 
references therein). 

The non-relativistic gravitational N-body problem, in particular, belongs to 
this class of systems. If one considers the motion of n + 1 particles (the Sun, 
with mass mo, and n planets with masses rrn, i = 1, . . . ,n) only affected by 
their mutual gravitational interaction, the corresponding equations of motion 
can be derived from the Hamiltonian 



where and Pi = rrii denote the position and momenta of the n + l bodies in 
a barycentric reference frame. Typically, however, the planets evolve around the 
central mass following almost Keplerian orbits, so that by an appropriate change 
of coordinates one can rewrite the Hamiltonian ([I]) as H = Hk + Hj , where in 
some sense \Hj\ <C l-ff/fl, or equivalently, as the sum of the Keplerian motion 
of each planet around the central mass and a small perturbation due to the 
gravitational interaction between planets. Jacobi and Heliocentric coordinates 
constitute paradigmatic examples of canonical set of coordinates possessing this 
feature. Thus, the Hamiltonian written as H = Hk + Hj is a particular 
example of a near-integrable Hamiltonian system, i.e, it can be expressed as 



into account this special structure when designing integration methods to ap- 
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proximate its dynamics. The idea consists in constructing splitting schemes as 
compositions of the flows corresponding to H^ a \q,p) and H^(q,p), assuming 
that they are explicitly computable or sufficiently well approximated 0GI]- In 
fact, since the parameter e is small, it is possible to design methods which be- 
have in practice as high order integrators with less severe restrictions concerning 
the order conditions than the usual split into kinetic and potential energy. This 
approach was systematically pursued by McLachlan [TS], obtaining families of 
splitting schemes of order 2 and 4 which eliminate the most relevant error terms 
in s, and further analyzed by Laskar & Robutel |1(J| in the context of planetary 
motion. 

By incorporating the idea of processing, even more efficient schemes can be 
constructed for the Hamiltonian ^ [2j. In that case, both the kernel and the 
processor are taken as compositions of the flows associated with ijM and , 
so that the exactly symplectic character of the integration scheme is ensured. 
With this approach, all terms of first order in e in the truncation error expansion 
can be annihilated with the processor |13|. [25] . 

Although the symplectic methods developed in |10| and |12] for near-integra- 
ble Hamiltonian systems have proved their usefulness in long time integrations 
of the Solar System [9] , the design of new and more efficient higher order inte- 
grators would certainly be of interest for numerical simulations of its evolution 
over a larger time span, either by speeding up the algorithms or by providing 
more accuracy in the position of the different objects. Relevant examples where 
the new integrators could be useful include the numerical integration of the So- 
lar System for 60 million years backward in time to cover the Palaeogene period 
to determine insolation quantities of the Earth and calibrate paleoclimatic data, 
studies of the planetary orbits over several billion years, etc. [TT| . To this pur- 
pose, it is essential that the numerical solutions obtained are not contaminated 
by error accumulations along the integration and that the computations are 
done in a reasonable time. 

The purpose of this work is to present new families of symplectic splitting 
methods specifically designed for Hamiltonian systems of the form ^ appear- 
ing in many problems of dynamical astronomy, when one is interested in highly 
accurate results over a large time span. The schemes we propose will be useful 
in particular in the the long time integration of the Solar System, both in Ja- 
cobi and Poincare Heliocentric coordinates, and are more efficient than schemes 
designed in [TO] and [12] . Although they involve the computation of the ap- 
plication of more elementary flows per step than other methods, their smaller 
error terms allow to use larger steps, which results in more efficient schemes. 
Obtaining the new methods requires deriving previously the necessary and suf- 
ficient order conditions to be satisfied by the coefficients (which is done here in 
a systematic way) and then solving these polynomial equations to get the best 
solutions according with some appropriately chosen optimization criteria. This 
is worked out in sections [2] [3] and the appendix, whereas in section [4] we consider 
the application of the new schemes to the integration of the Solar System. The 
new methods obtained in section [3] are suitable to be applied when using Jacobi 
coordinates, but not very appropriate for the integration of the Solar System in 
Poincare Heliocentric coordinates. In subsection 4.2 we show that, in the later 
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case, it is more convenient to use splitting methods whose coefficients satisfy an 
additional polynomial equation, and we construct several schemes of this class 
which are as efficient as those constructed in section [31 

It is worth stressing that, while the main motivation of this work is the long 
time integration of Hamiltonian problems arising in dynamical astronomy, and 
in particular in planetary systems, the new symplectic splitting methods ob- 
tained here can also be applied to more general perturbed differential equations 
arising in different fields when high accuracies are required. 



2 Order conditions 
2.1 Preliminaries 

To establish the framework for the construction and analysis of the new families 
of integrators, we consider a generic differential equation of the form 

a/ = fW(x) + ef®(x), x(0) = x eR D , (3) 

where \e\ <C 1 and each part 

x' = f [a] {x), x' = ef lb] (x) (4) 

is exactly solvable (or can be numerically solved up to round off accuracy) with 
solutions 

x(t) = ^(x ), x(t) = ^(x ) 

respectively, at t = r, the time step. If we denote by (p T (xo) the exact solution of 
it is well known that ijj T = ip$ ° ^fr^ provides a first-order approximation, 
i.e., i/j t (xo) = ip T (xo) + 0(t 2 ) and that higher order approximations can be 
obtained by taking more compositions in ip T , 

for appropriately chosen coefficients aj,fej. The splitting method ip T is said to 
be of order r if for all x G M. D , 

4> T {x) = ifrix) + 0(r r+1 ) as r -»• 0. (6) 

It is straightforward to check that the method is at least of order 1 for arbi- 
trary problems of the form ^ if and only if the coefficients o« , b{ satisfy the 
consistency condition 

s+l s 
1=1 i=l 

We are mainly interested in symmetric methods, that is, integrators verify- 
ing V'-t = ipr i or equivalently a s+ 2-i = ai, = 6j (so that the composi- 
tion ([5]) is left-right palindromic). In that case, they are automatically of even 
order. In particular, if a symmetric method satisfies the consistency condition 
([7]), then it is at least of order 2. 



4 



Since the last flow (pa} +1T can be concatenated with the first ipair at the 
next step when scheme ^ is iterated, the number of flows (pfi and <p$ per 
step is precisely s. This number is usually referred to as the number of stages 
in the composition. 



2.2 Deriving the order conditions via the BCH formula 

The conditions that the coefficients a,i,bi must satisfy for a splitting method 
to be of order r (the so-called order conditions) can be conveniently derived 
by considering series of linear differential operators. We denote by A and B 
the Lie operators associated with and respectively. For each smooth 
function g : MP — > M, Ag and B g are smooth functions defined as 



Ag{x)=± 



g{^{x)), Bg{x)=^ 



T=0 



for each x S M. D , that is, 

Ag(x) = fW(x)-Vg(x), B g(x) = f®(x) ■ Vg(x). (8) 

The near-integrable Hamiltonian system ^ corresponds in this general frame- 
work to considering equation ^ with 

x = (q,p), f [a] (x) = JVH^(q,p), and f [b] (x) = JVH®(q,p), 

being J the canonical symplectic matrix. Therefore, for each smooth function 
g one has 

Ag = f^-Vg = T—^--—^- 

j d Pj 'hi ih u <h>j ' 

and a similar expression for B g. 

It is well known that for any smooth function g, the r-flow of ^ satisfies 

g^ T (x)) = ^ A+ ^g{x), 

where e T ( A + eB ) [ s defined as a series of linear differential operators 



OO L, 



e T{A + sB) = Y^L- {A + eB) k . 

The same is true for each part in ([3]): 

g(^(x)) = e^ A g(x), g(^ ] (x)) = e^ B g(x). (9) 

Analogously, for the integrator i)j t in ([5]), one has 

g(ip T (x)) = V(r)g(x), 
where ^(t) is a series of linear differential operator defined as 

*&(r) = e airA e blTeB ■ ■ ■ e aaTA e baTeB e a s+i Tj4 _ (\ q\ 
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Notice that the exponentials of Lie derivatives in (10) appear in the reverse 
order with respect to the maps in the integrator ([5]). 

One of the standard ways of deriving the order conditions for splitting meth- 
ods is the following. By applying repeatedly the Baker-Campbell-Hausdorff 



(BCH) formula [22] to the factorization (10) corresponding to a consistent (i.e., 
satisfying ([7])) splitting method, one is able to express ^(r) as the formal ex- 
ponential of only one operator: 



where 



^( r ) _ e T(A+eB+E(r)) ^ 

E(t) = TE Pab [A,B]+T 2 ep aba [[A,B],A]+T 2 e 2 Pabb [[A,B},B]) 

+r 3 e Pabaa [[[A, B],A],A] + r 3 e 2 Pabba [[[A, B], B},A] (11) 
+T 3 e 3 Pabbb [[[A,B},B},B} + 0(r 4 ). 

Here the symbol [A, B] stands for the commutator of the Lie operators A and 
B, and p a b,Pabb,Paba,Pabbb, ■ ■ ■ are polynomials in the parameters a^bi of the 
splitting scheme. In particular, 

1 A, lA, , ,1 



.-, X^iCi, Paba = jjJ^W 1 12 
i=l i=l 



where 



1,2,..., s (12) 



and c s +i = 1. The integrator is of order r if E(t) in (11) is of size 0(r r ), 
so that ^(r) agrees with the series of linear operators e T ( A + eB ) Q f the exact 
flow up to terms of size 0(t v ). In consequence, the order conditions read 
Pab = Pabb = Paba = ' ' ' = up to the order considered. For symmetric methods, 
E(t) only involves even powers of r, that is, p w = for any word w with an 
even number of letters in the alphabet {a, b}. 



In (11) we have considered the classical Hall basis associated to the Hall 



words a, b, ab, abb, aba, abbb, abba, abaa, . . . [T7]. The coefficients p w in ( 11 ) cor- 
responding to each Hall word w can be systematically obtained using the results 
in |16j in terms of rooted trees and iterated integrals. An efficient algorithm 



(based on the results in [TE]) of the BCH formula and related calculations that 



allows one to obtain expression (11) up to terms of arbitrarily high degree is 
presented in [3]. 



2.3 Generalized order 

We are particularly interested in the manner in which the local error ip T {x) — 
(f T (x) decreases as e — > 0. For instance, from the results in the precedent 
subsection, it is clear that for any consistent symmetric method the local error 
satisfies ip T (x) = (p T (x) + C(er 3 ). Alternatively, 

*(r) -e TiA+£B) = 0{et 3 ) as (r,e) (0,0). 
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If in addition p a bb = in ( [TT] ) , then 



i> T {x) = <p T (x) + 0(er 5 + e 2 r 3 ) as (r, e) (0, 0) . 

In that case, we say that such a method is of (generalized) order (4,2). More 
generally, we will say [H] that a splitting method is of generalized order (r±, r2, ■ ■ ■ , r m ) 
(where r\ > r2 > • • • > r m ) if 

i> T {x) = ip T {x) + 0(er ri+1 + eV 2+1 + • • • + e m r rm+1 ). 



2.4 Generalized order conditions 

The conditions that the coefficients a% , hi must satisfy for a splitting method to 
be of a prescribed (generalized) order (n, r2, • ■ . , r m ) can be obtained, of course, 
by computing the polynomials p w in expression ( |11[ ) with the BCH formula and 
then equating each term to zero up to the considered order. Thus, in particular, 
a consistent symmetric scheme of order (6, 2) requires that p a f, a = p a baaa = 0. 

There exist, however, other more systematic procedures to derive these order 
conditions. In what follows, we present a strategy that allows us to get in a 
direct way a set of necessary and sufficient independent order conditions for 
generic splitting methods. 

As a first step, we consider Z(t) = e T ( A + £B ) e ~ TA , which is the formal 
solution of the initial value problem 

^Z(t)=sZ(t)C(t), Z(0) = I, (13) 

where 

oo 

C( T ) = e TA Be~ TA = r n ~ x C n , (14) 

n=l 

with 

d = B, C n = T - i —-\A,[A,...,[A,B})}, n>l. 
(n - iy. v v / 

n— 1 times 

On the other hand, applying repeatedly the identity e TA e hB e~ rA = e hC ^ to 
eq. ([To]) and taking into account (12), we arrive at 

*(r) = Z(r)e Cs+1TA , where Z(t) = e £hlTC{ciT) ■ ■ ■ b ^ s tC{c s t)_ ^ 

Notice that, if the splitting method is consistent, then c s +i = 1. We thus have 
that a splitting method is of order (n, r%, . . . , r m ) if and only if c s+ \ = 1 and 

Z(t) - Z(t) = 0(er ri+1 + eV 2+1 + • • • + e m r r ™ +1 ). (16) 

We then expand both Z(t) and Z(t) as power series of e and compare their 
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coefficients. First, applying Neumann iteration to (13) we get 
Z(t)-I = e I Z( Sl )C{si)dsi 



J o 

= e r c( S1 )d S1 +£ 2 r r c{ S2 )c{s 1 )ds 1 ds 2 

Jo Jo Jo 

+ £ 3 [ I 1 f 2 C(s 3 ) C{s 2 ) C(si) d Sl ds 2 ds 3 + --- 
Jo Jo Jo 

= g £ * . £>, (ft + ■■■ + *)•••«.+ « A ° h ■ ■ ■ Ch ' 

where in the last equality we have introduced explicitly the expression for C (r) 
given by (14). On the other hand, by expanding the exponentials of Z{r) in 
(15), we have 

s 

Z(t)-I = T£Y J hC(c l T) 
1=1 

+T 2 e 2 (^|c( Ci r) 2 + X; E hbjC^Cicjr) ] +••• 
y«=i i=i i=i+i / 

= £tV E ^^^r)--^) 

' * ' * CTr 



k>l l<il<-<ifc<s 



where 



E^ fe E jk ( E ^^i^-^r 1 )^---^ 

fc>i ji,-.j'fe>i \i<n<-<*fe<« %1 "" lk J 



1 if ii < ••• < i k , 



1 ., . . . . . 

a h-ik = J a it+i-ik 11 z i = ' ' ' = H < U+l S ■ ■ ■ < ^k■ 

In this way a splitting method is of order (n, . . . , r m ) if and only if 



Kil < L < it < s ^r-ik n ' 0l + -+i*)-(7l+Ja)jl 



for each /c = 1, ... ,m and each multi-index (i.e., fc-tuple of positive integers) 

(j'l, ...,jk) such that ji H \-jk<r k - 

Conditions (17) (one condition for each multi-index) have been obtained 
in [21] in the context of order conditions of splitting operators for unbounded 
operators A and B. Nevertheless, such order conditions are not all independent. 
For instance, it can be checked that if condition ( |17[ ) holds for the multi-indices 
(1,2), (2), and (1), then the condition for (2,1) is also fulfilled. That kind 
of dependencies are a consequence of the fact that both Z(r) and Z[t) are 
exponentials of Lie series in the non-commuting indeterminates C\ , C 2 , ■ ■ ■ . A 
set of independent order conditions can be obtained (by virtue of Theorems 3.2 



S 



and 6.1 in [T7]) by considering a particular subset of multi-indices, the so-called 
Lyndon multi-indices. Let us consider the lexicographical order < (i.e., the 
order used when ordering words in the dictionary) on the set of multi-indices. 
A multi-index (ii, . . . , i m ) is a Lyndon multi-index if (ii, . . . , i^j < {ik+i-, ■ ■ ■ , im) 
for each 1 < k < m. For instance, the subset of Lyndon multi-indices {j\ , . . . , jk) 
such that ji + • • • + jk < 5 is 

{(1), (2), (3), (4), (5), (1,2), (1,3), (1,4), (2, 3), (1,1,2), (1,1,3), (1,2,2), (1,1, 1,2)}. 

Taking into account these considerations, we finally arrive at the following 
result. 

Theorem 1 A splitting method of the form |5|) is of generalized order (n, . . . , r m ) 
if and only if c s +i = 1 and (17[ ) holds for k = 1, . . . , m and each Lyndon multi- 
index (ji, . . . ,jk) such that ji + • • • + jk < rk- For symmetric methods, only 
Lyndon multi-indices (ji, . . . , jk) with odd ji + • • • + jk need to be considered. 



For illustration, in Table [T] we collect explicitly conditions (17) corresponding 
to some particular multi-indices, whereas in Table [2] we specify which particu- 
lar Lyndon multi-indices one has to consider, or equivalently which conditions 



(17) must hold for each consistent symmetric splitting method of the given 



generalized order, according to Theorem [T] 



Multi-index 


Condition 


U), i>i 


i=i J 

Yl 2% Ci + Yl bib i c i = 3 

i=l Ki<j<s 


(1,2) 


(1,4) 


S 1 1 

i=l l<!<j<S 


(2,3) 


s 1 1 

E h J>M = w 

i=l l<«<i<s 



Table 1: Generalized order condition associated with each Lyndon multi- index. 

To this point some remarks must be done. The set of order conditions 
given by Theorem [T] is completely equivalent to the order conditions that can 



be obtained by following the standard approach described in Subsection 2.2 
On the one hand, the free Lie algebra C(C\, C2, C3, . . .) generated by the non- 
commuting indeterminates Ci, C2, C3, . . ., admits a basis (the Lyndon basis [TT] ) 
in one-to-one correspondence with the set of Lyndon multi-indices. Clearly, if 
instead of directly comparing the series expansions of Z(r) and Z{t) as above, 
we compare the formal logarithms log(Z(r)) and log(Z(r)), we could obtain 
one order condition per element in the Lyndon basis. On the other hand, the 



approach in subsection 2.2 gives one order condition per element in a basis of 
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Generalized order 


Lyndon multi-indices 


(2n,2) 


(3),(5),...,(2n-l) 


(8,4) 


(3), (5), (7), (1,2) 


(10,4) 


(3), (5), (7), (9), (1,2) 


(8,6,4) 


(3), (5), (7), (1,2), (1,4), (2,3) 


(10,6,4) 


(3), (5), (7), (9), (1,2), (1,4), (2,3) 



Table 2: Lyndon multi-indices corresponding to consistent symmetric splitting meth- 
ods of a given generalized order. 

the free Lie algebra C(A, B) generated by the noncommuting indeterminates 
A and B, and Lazard elimination theorem [T7] shows that, as vector spaces, 
the direct sum of £(Ci, C2, C3, . . .) with the linear span of A is isomorphic to 
C{A, B) (sending C x , C 2 , C 3 , . . . to B, [A, B], \[A, [A, B]},... respectively). 

3 New numerical schemes 

There are two different types of symmetric composition schemes ([5]): one in 
which the first and last flows correspond to the A part (and thus appropriately 
called ABA composition), 

ABA: <pjg T o $1 o o • • • o <pW. o v W t o V \< T (18) 

and the other in which the role of tpf and ipip is interchanged (BAB composi- 
tion): 

BAB: o o v fl o • • • o ^ o V W. o ^ (19) 

Notice that both types of composition are closely related: an s-stage BAB 
method is just an (s + l)-stage ABA scheme with a± = 0, so that to construct 
BAB methods one has to solve the same order conditions as for ABA com- 
positions. Although A and B are qualitatively different here, and therefore 
both types of composition may lead in principle to integrators with different 
performances, in practice, and for the examples analyzed, we have not found 
substantial differences, so that in what follows we only consider ABA methods 
for clarity in the presentation. 

Constructing particular methods requires solving polynomial equations (e.g. 
the order conditions of Table [T] for methods of Table [2]) , a problem whose com- 
plexity grows enormously with the number of equations and variables involved. 
This task can be handled by computer algebra systems when this number is 
relatively low. In that case one is able to get all the solutions and select the 
one that verifies some previously fixed optimization criterion, such as minimiz- 
ing error terms at higher orders and the the sum of the absolute value of the 
coefficients. In practice, we have followed this procedure when there are no 
free parameters and the number of equations to be solved is at most seven. 
When this number is larger than seven or there are additional parameters, an- 
other strategy based on homotopy continuation methods has been applied. In 
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the appendix we provide a detailed treatment of the procedure for a particular 
method. 



3.1 New methods in the ABA class 

Symmetric schemes of generalized order (2n, 2) can be obtained just by solving, 
in addition to consistency, the order conditions corresponding to the Lyndon 
multi-indices (3), (5), ... , (2n — 1) (first line in Table [I]). These equations result 
from approximating the integral J Q C(s)ds in the expression of Z{t) by the 
quadrature rule 

i=l i=l j>l 

in the expansion of Z(t). Equivalent order conditions were previously derived 
in \10\ [T2"l 113] . and so the same methods are obtained here. Methods in this 
family have all coefficient positive and good stability properties. In the tests 
carried out in this paper we will take the most efficient ABA scheme of order 
(8, 2) for comparison, which we denote by ABA82. 

Generalized order (10,4). According with Table [2j there are five order 
conditions in addition to consistency ([7]), for a total number of seven equations 
to be satisfied by the coefficients. In consequence, the minimum number of 
stages is six. A more efficient method can be obtained, however, by taking an 
additional stage and choosing the corresponding free parameter to reduce the 
error terms at a higher order. The sequence of coefficients is 

«i h a 2 b 2 a 3 b 3 a 4 6 4 a 4 b 3 a 3 b 2 a 2 h a\ (20) 

and their values are collected in Table [3] (method denoted by ABA104). Observe 
that, as expected, one of the a, and one of the bj coefficients are negative (it is 
know that this feature is unavoidable for any splitting method of order higher 
than two [5] 1191 [20]. but they have a relatively small absolute value. 



Generalized order (8,6,4). Now we have, in addition to consistency, six 
order conditions, for a total of eight equations, so that the minimum number 
of stages is seven. Hence the sequence of coefficients for the resulting methods 
is also (20). There are 30 real solutions, and the one referred to as method 
ABA864 in Table [3] minimizes the sum of the absolute values of the coefficients. 



Generalized order (10,6,4). An additional stage is required in this case to 
verify the order condition associated with multi-index (9) in Table [T] Therefore, 
the minimum number of stages is eight, with sequence 

oi h a 2 b 2 a 3 b 3 a 4 6 4 a 5 6 4 a 4 63 a 3 b 2 a 2 b\a\. (21) 

By following the construction strategy exposed in the appendix, we have ob- 
tained several solutions for the order conditions. Among those possessing rea- 
sonably small coefficients, we have selected the solution with the smallest lead- 
ing terms of the local error. This corresponds to method ABA1064 in Table [3} 
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id order stages di,bi 









<Xl 


= 





. 047067100645972506129478876372 








«2 


= 


0, 


. 184756935417088106924737619370 








0,3 


— 


0. 


. 282706005679836205324361656554 


ABA 104 


(10,4) 


7 


«4 


— 


-0. 


.014530041742896818378578152296 


bi 


= 


0. 


. 118881917368197019945350395085 










= 


0. 


. 241050460551501565744166786590 








I 

03 


= 


-0 . 


. 273286666705323806054311398166 








I 

04 




. 


. 826708577571250440729588432981 








(11 


= 


0, 


.071133426498223117777938730006 








ai 


= 


0. 


. 241153427956640098736487795326 








0,3 


= 


0. 


.521411761772814789212136078067 


ABA864 


(8,6,4) 


7 


"4 


= 


-0. 


. 333698616227678005726562603400 


bi 




0. 


183083687472197221961703757166 








l>2 




0. 


.310782859898574869507522291054 








b 3 




-0. 


.026564618511958800697212137916 








bi 




0, 


. 065396142282373418455972179391 








Ol 




0, 


. 038094497422412195456975322308 








«2 




0. 


. 145298716116913749294020072660 








03 




0. 


. 207627695725541250716205611324 








"4 




0. 


.435909703651526159223154862401 


ABA1064 


(10,6,4) 


8 


as 




-0. 


. 653861225832786709380711737390 








bi 




0. 


. 095858880837075210610771503771 








&2 




0. 


. 204446153142998780680507783916 








6 3 




0, 


.217070347978991101714338592430 








In 




-0. 


017375381959065093005617880118 



Table 3: Coefficients for ABA symmetric splitting methods of generalized order 
(10,4), (8,6,4) and (10,6,4). 



3.2 A simple example 

To illustrate the efficiency of these schemes we take the perturbed Kepler prob- 
lem with Hamiltonian 

where r = \/qf + q%- This Hamiltonian is a first approximation used to de- 
scribe the dynamics of a satellite moving into the gravitational field produced 
by a slightly oblate spherical planet and whose motion takes place in a plane 
containing the symmetry axis of the planet |15| . 

We consider this simple problem to test the relative performance of the 
methods obtained in this work in comparison with schemes presented in |10|I12|. 
The following schemes are used: 

• ABA82: The 4-stage (8,2) ABA method given in [10} H~2]. 

• ABA84: The 5-stage (8,4) ABA method of [12]. 

• ABA104: The 7-stage (10,4) method given in Table [3} 

• ABA864: The 7-stage (8,6,4) method of Table [3} 
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• ABA1064: The 8-stage (10,6,4) method whose coefficients are collected in 
Table! 

For the numerical experiments we take as initial conditions qi = 1 — e, 
Q2 = 0, p\ = 0, pi = + e)/(l — e), with e = 1/4, which would correspond 
to the eccentricity of the unperturbed Kepler problem. For this system, the 
strength of the perturbation depends both on the choice of the small parameter, 
e, and the initial conditions. We integrate along the interval t £ [0, 10000] 
and compute the averaged error in energy as well as the averaged error in 
position and momenta (measured in the 2-norm) of the numerical solutions 
evaluated at = 20 • k, k = 1,2,..., 500. We take as the exact solution an 
accurate approximation obtained using a high order method with a sufficiently 
small time step. This numerical test is repeated several times for each method 
using different time steps (changing the computational cost for the numerical 
integration). Finally, we plot the average errors versus the time step scaled by 
the number of stages per step, i.e. t/s, in double logarithmic scale, to show the 
error versus the computational cost (the cost is inversely proportional to t/s, 
and the best methods should provide a given accuracy with the largest value of 
t/s). 

Figure [l] shows the results obtained for e = 10 2 ,10 3 . In the left panel 
we show the average error in positions and momenta, and in the right panel 
we repeat the same experiments, but in this case we measure the average error 
in energy. Notice how the new methods collected in Table [3] are clearly more 
efficient than ABA82 and ABA84. 



4 Application to the integration of the Solar System 

In this section we illustrate how the new families of methods proposed here 
behave when they are used in the numerical integration of the simplest model 
of the Solar System, i.e., a main massive body (the Sun) and a set of particles 
(the planets) orbiting the Sun following almost Keplerian trajectories. It is not 
our intention to carry out a comprehensive treatment of this problem, but rather 
to check the performance of the new methods and compare them with other 
well established schemes designed for near-integrable Hamiltonian systems such 
as those presented in [10] and |12j . We instead refer the reader to the reference 
[3], where this issue is handled in much more detail. It is worth remarking 
that all the simulations in this section have been done using and extended real 
arithmetics and that we use compensated summation during the evaluation of 
the internal stages of the schemes. 

As stated in the Introduction, integrating numerically the gravitational N- 
body problem requires first to choose a convenient set of canonical coordinates. 
Two widely used coordinate systems where the corresponding Hamiltonian ([!]) 
adopts the form ([2]), amenable to the application of the integration schemes 
developed in this work, are Jacobi and Heliocentric coordinates. In the former, 
the position of each planet is taken relative to the barycenter of the previous i 
bodies, whereas in the later the position of the central star (the Sun) is chosen 
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-2 -1.8 -1.6 -1.4 -1.2 -1 -2 -1.8 -1.6 -1.4 -1.2 -1 



Figure 1: Averaged error in positions and momenta (left panels) and averaged 
error in energy (right panels) versus the scaled time step, t/s, in a double 
logarithmic scale for the numerical integration of the Hamiltonian system (22) 
along the time interval t G [0, 10000] and measured at times = 20 ■ k, k = 
1,2, ...,500. 



as the center of coordinates and the position of each planet is taken with respect 
to the Sun. 



4.1 Integration in Jacobi coordinates 

The positions are denoted by Vq, Vj, . . . v n; where vo is the center of mass of 
the system. In terms of the original coordinates qi of they read [8] 

vo = (m-oqoH \-m n q n )/r] n (23) 

i-l 

Vj = qi - (^mjcij)/^-!, j = i,...,n (24) 
whereas the corresponding canonical momenta are given by 

i-l 

vo = Po H HPn, Vj = (ryj_ipi - m, ^qj)/^. (25) 

j=0 
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Here rji 
form 



rrin. 



In this set of coordinates the Hamiltonian 




takes the 



(26) 



where A« stands for the distance between bodies i and j and therefore can be 
expressed in terms of Vj and Vj. In consequence, it belongs to the class Q 
with the perturbation depending only on positions and therefore all the 
integration methods presented in the previous section can be directly applied. 



Mer - Ven - Ear - Mar (Jacobi Coord) 



Jup - Sat - Ura - Nop (Jacobi Coord) [ short ] 



Mercury to Neptune (Jacobi Coord) 




ABA82 ■ 
ABA84 
ABA104 ■ 
ABA864 ■ 
ABA1 064 




ABA82 
ABA84 
ABA104 ■ 
ABA86 
ABA1064 




Figure 2: Comparison between the ABA splitting schemes (10,4), (8,6,4) and (10,6,4) 
of Table [3] against ABA82 and ABA84. From left to right: the four inner planets, the 
four outer planets and the whole Solar System. The x-axis represents the (inverse of 
the) cost r/s, and the y-axis the maximum energy variation for one integration with 
constant step size r. Both are in logarithmic scale. 

In Figure [2] we can see the results on the N-body problem in Jacobi coordi- 
nates for different planetary configurations. From left to right we have the four 
inner planets (Mercury, Venus, Earth and Mars), the four outer planets (Jupiter, 
Saturn, Uranus and Neptune) and the eight planets in the Solar System (Mer- 
cury to Neptune). The initial conditions and mass parameters have been taken 
from the JPL Solar System ephemerides DEA405 (http : //ssd. jpl .nasa. gov/). 

We have integrated the same initial conditions with the different ABA split- 
ting schemes presented in section [3j using different step sizes r. For each step 
size (Tj = 1/2* for i = 1, 15) we have computed the trajectory over niter = 
10 5 evaluations of the integration scheme (i.e. if r = 0.5, then the final time is 
tf = 50000years). For each trajectory we plot the maximum variation in energy 
along the trajectory versus the inverse of the computational cost, t/s, both in 
logarithmic scale. Notice that we do not integrate the each initial condition to 
the same final time but this should not affect to the results as there is no drift 
on the energy during the computation. 
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As we can see for the four inner planets and the whole Solar System (Figure [2] 
left and right respectively), the performance of the different ABA schemes in 
Table [3] is equivalent to ABA84 constructed in |12j and much better than method 
ABA82 given in |10} I12j. However, if we focus on the simulations done for the 
four outer planets (Figure [2] middle) we see that the ABA schemes of orders 
(8,6,4) and (10,6,4) are much better than ABA84. 

Despite the fact that the size of the perturbation for the different planetary 
configurations presented can be very different [1] , we believe that the difference 
between the performance of the schemes for the inner and the outer planets in 
the Solar System is mainly due to Mercury. Its fast orbital period and relatively 
high eccentricity are the main limiting factors when one tries to improve the 
efficiency of the higher order schemes. Notice that for the inner planets the 
orbital period is much shorter than for the outer planets, and so, according 
with [23j , this imposes a restriction on the step size to be used and the number 
of evaluations per orbital period. 

4.2 Integration in Poincare Heliocentric coordinates 

In this set the coordinates rj are the relative positions of each planet with 
respect to the Sun: 

r = q , rj = qj-q , i = l,...,n (27) 

whereas the conjugate momenta read 

?0 = Po H \-Pn, f j = pi (28) 

and the Hamiltonian ([!]) is given by [8] 

gHe = f 1 n- .112 mo + m c momj \ | c m i m A ? (29) 

i=l " " ' 0<i<j<n J 

where Ajj = ||rj — r^H for i, j > 0. 

Heliocentric coordinates have the advantage, in comparison with Jacobi co- 
ordinates, that adding a new body to the model does not change the origin 
and the new Hamiltonian is easily updated. On the other hand, the pertur- 
bation depends on both positions and momenta and is not integrable by 
itself. Hence, when considering the splitting ([3]), the equation x 1 = f^(x) is 
not exactly solvable and the previous setting has to be generalized to cope with 
this situation. However, we observe that is the sum of two terms, one 
depending only on positions and the other depending only on momenta 

so each equation x' = /™(ar), with p bi \x) = JV (x) , is exactly solvable 
with flow ifr^ ■ One possibility consists in taking splitting methods of the form 

W = o *1 o rfj o ^H T o • • • o #1 o $1 o <pM T 
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and obtaining the appropriate coefficients Oj, bi, Ci satisfying the required order 
conditions. These can be derived, for instance, by analyzing the free Lie algebra 
generated by the three Lie derivatives corresponding to each piece of the Hamil- 
tonian. The number of order conditions grows rapidly with the order, however, 
in comparison with splitting schemes involving only two parts. In particular, 
(10,6,4) methods require, in addition to consistency, 22 order conditions. 

A much simpler alternative consists in replacing the exact flow tp^} T in the 

splitting method (5) by the second order symmetric approximation tj$} pro- 
vided by the leapfrog scheme, 

~V>] [h] [&a] [61] /onA 

since cpf] = (pf: ' +0(r 3 e 3 ). Then, the corresponding series of linear differential 
operators associated with this composition reads 

This implies that the additional error term introduced in the composition (pj) by 
the approximation (30) can be made to be of order 0(r 5 e 5 ) just by considering 



the additional order condition 

s 

£> 3 = 0. (31) 

i=l 

In this way, symplectic splitting schemes for Heliocentric coordinates can be 
designed along the same ways as in section 3 with an additional stage to verify 



equation (31). A consistent (10,6,4) method of this type has to satisfy only 8 



order conditions (instead of 22). It should be stressed that methods satisfying 



the additional restriction (31) are valid for general problems of the form 

x > = fW(x) + efM(x) + ef®(x), 
and not only for planetary problems in Heliocentric coordinates. 



Notice that ABA-type compositions (18) are clearly more convenient now, 
since the last stage of the method in the current step can be concatenated with 
the first stage at the next step. This is not possible with BAB compositions 



(19), because 

By following a similar strategy as for the methods collected in section |3j 
we have constructed symmetric schemes within this family of generalized order 
(8,4), (8,6,4) and (10,6,4), all of them involving the minimum number of stages. 
The corresponding coefficients are collected in Table |4j For schemes (8,4) we 
have found all the real solutions and selected the solution that minimizes the 
sum of the absolute values of the coefficients (method ABAH844). The procedure 
for constructing method ABAH1064 is detailed in the appendix, and a similar 
strategy has been used to build scheme ABAH864. 
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id 


order 


stages 








7 

ai, bi 








CI i 




. 


. 27414026894340187616405654402 








(1 2 




-0 


1 075684384401 642306251 1 052968 








(I Q 




-0 . 


. 04801850259060169269119541721 


ABAH844 


(8 4) 


6 


Old 




. 


. 76289334417472809430449880574 








bi 




. 


. 64088579516251271773224911649 








u 2 




-0 . 


. 85857544895678285658812832469 








bs 


= 


0. 


. 71768965379427013885587920820 








(I I 

1 




. 


. 06810235651658372084723976682 








Cl-i 




. 


. 251 13603872210332330728295804 












-0 . 


. 07507264957216562516006821767 








(7. 4 




-0 . 


. 00954471970174500781148821895 


ABAH864 


(8,6,4) 


8 


as 


= 


0. 


. 53075794807044717763406742353 








bi 




. 


. 16844325936189545343103826977 








u 2. 




. 


. 42431771737426772243003516574 








bi 




-0 . 


. 58581096946817568123090153554 








bi 


= 


0. 


. 49304999273201250536982810002 








ai 




. 


. 047319086976533822704043 1 1796 








a 2 




0. 


. 26511052357487851595394800361 








C13 




-0. 


. 00997652288381 124084326746816 








04 




-0. 


. 05992919973494155126395247987 


ABAH1064 


(10,6,4) 


9 


tl5 
6! 




0. 
0. 


. 25747611206734045344922822646 
. 11968846245853220353128642974 








b 2 




0. 


. 37529558553793742504201285376 








6, 




-0. 


. 46845934183259937836508204098 








h 




0. 
0. 


. 33513973427558970103930989429 
. 27667111912108009750494572633 



Table 4: Coefficients for ABA symmetric splitting methods of generalized order (8,4), 
(8,6,4) and (10,6,4) especially adapted to be used with Heliocentric coordinates. 



In Figure [3] we can see the results on the N-body problem in Heliocentric 
coordinates for different planetary configurations. As for Jacobi coordinates, 
from left to right we have the four inner planets (Mercury, Venus, Earth and 
Mars), the four outer planets (Jupiter, Saturn, Uranus and Neptune) and the 
eight planets in the Solar System (Mercury to Neptune) . 

As before, we have integrated the same initial conditions with the different 
ABA splitting schemes presented in Table |4j using different step size r. For 
each step size (r^ = 1/2* for i = 1, . . . , 15) we have computed the numerical 
trajectory over niter = 10 5 evaluations of the scheme, and for each trajectory 
we plot the maximum variation in energy along the trajectory versus the inverse 
of the computational cost, t/s, both in logarithmic scale. 

Notice that in the case of the four inner planets (Figure [3] left) the perfor- 
mance for the different ABA schemes in Table |4] are better than method ABA82 
given in |XQ^ I12j . but there is not much different between the ABA schemes 
ABAH844, ABAH864 and ABAH1064, especially designed for Heliocentric coordi- 
nates. Nevertheless, if we look at the results for the four outer planets and the 
whole Solar System (Figure [3] middle and right respectively) we can see that 
ABAH864 and ABAH1064 show better results than ABAH844. 

Finally, it is worth mentioning that with the ABA schemes of orders (8,6,4) 
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and (10,6,4) there is no longer a big difference between Jacobi and Heliocentric 
coordinates. Looking at the results for the two sets of coordinates systems 
for the whole Solar System (Figures [| and [3] right) we notice a big difference 
between the optimal step size if we take the ABA (8,4) scheme. Where the 
optimal cost is ~ 1CP 4 in Heliocentric coordinates against ~ 10~ 3 for Jacobi 
coordinates. While for the ABA (8,6,4) and ABA (10,6,4) the optimal cost for 
both set of coordinates is ~ 10~ 3 . This is important as Heliocentric coordinates 
are easier to handle. 



Mer - Ven - Ear - Mar (Helio Coord) Jup - Sat - Ura - Nep (Helio Coord) Mercury to Neptune (Helio Coord) 




-5 -4 -3 -2 -1 -5 -4 -3 -2 -1 -5 -4 -3 -2 -1 



Figure 3: Comparison between ABA schemes of order (8,4), (8,6,4) and (10,6,4) 
especially designed for Heliocentric coordinates and ABA82. From left to right: 
the 4 inner planets, the 4 outer planets and the whole solar system. The x-axis 
represents the (inverse of the) cost r/s, and the y-axis the maximum energy 
variation for one integration with constant step size r. Both are in logarithmic 
scale. 



5 Concluding remarks 

In reference [10] , symplectic splitting methods of generalized order (2n, 2) up 
to n = 5 (first described in [12]) were systematically derived and tested on the 
Sun- Jupiter-Saturn system over 25000 years in Jacobi coordinates, observing 
an improvement in the accuracy with respect to the leapfrog integrator by 
several orders of magnitude at the same computational cost. Methods in this 
family have all the coefficients positive and good stability properties. Scheme 
(8, 2) in particular has been used in several long term simulations of the whole 
Solar System (e.g., [12 [H]) and corresponds to method ABA82 in the examples 
reported here. All the tests carried out in |10j showed that the error term 
e 2 r 3 was the main limiting factor in the performance of the integrators, so 
the natural question was whether schemes of higher order (and thus already 
involving some negative coefficients) could be useful for integrating planetary 
N-body problems. 

As a matter of fact, methods of order (8, 4) obtained in |12] do improve 
the performance of ABA82 for this problem in Jacobi coordinates, as the ex- 
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periments reported here and in [4J show. In this work we have pursued this 
line of research and constructed new families of higher order splitting methods 
specifically oriented to the numerical integration of near-integrable Hamiltonian 
systems, and in particular for planetary N-body problems, both in Jacobi and 
Poincare Heliocentric coordinates. To do that, first we have derived explicitly 
the set of independent necessary and sufficient order conditions that splitting 
methods must verify to achieve a certain order of accuracy and then we have 
solved these equations. A non-trivial task that requires the use of homotopy 
continuation techniques and optimization criteria to select the most appropriate 
solution. 

Although the new methods involve some negative coefficients, and thus one 
could think that their numerical stability might be compromised, they have been 
selected to minimize the error terms at higher orders and the sum of the absolute 
values of the coefficients. As a result, the size of the negative coefficients of our 
new methods is relatively small. In any case, the experiments reported here 
clearly indicate that the new methods of order (8,6,4) and (10,6,4) achieve 
accuracy up to round off error with larger step sizes than 2nd-order schemes. 

For N-body planetary systems, we have generalized these methods to be 
used also in Poincare Heliocentric coordinates just by adding an extra stage 
to the original composition designed for Jacobi coordinates. The new methods 
thus constructed are in fact valid for generic differential equations of the form 
x > = /H( x ) + efW(x) + e/M(x), where |e| <C 1. 

Our simulations show that the efficiency of the new integrators is essentially 
similar in both types of coordinates. We believe this result is worth remarking, 
since canonical Heliocentric variables provide very often a more convenient for- 
mulation of the problem. The improvement of the new integrators presented 
here (in particular, methods of order (8, 6, 4) and (10, 6, 4)) with respect to pre- 
vious schemes is most notably exhibited when they are applied for the numerical 
integration of the outer planets. When the whole Solar System is considered, 
although methods of order (8,6,4) and (10,6,4) still provide the best results, 
it is Mercury with its relatively high eccentricity and fast orbital period which 
constitutes the main limiting factor in all simulations. 

When designing the new methods of generalized order (8, 6, 4) and (10, 6, 4), 
both in Jacobi and Heliocentric variables, we have only considered compositions 
with the minimum number of stages to solve all the order conditions. It might 
be the case, as in other contexts, that introducing more stages with additional 
free parameters could lead to more efficient schemes. We intend to explore this 
possibility and eventually collect the new methods obtained, both in the ABA 
and BAB classes, in our website. 

Acknowledgements 

The work of SB, FC, JM and AM has been partially supported by Ministerio de 
Ciencia e Innovacion (Spain) under project MTM2010-18246-C03 (co-financed 
by FEDER Funds of the European Union), whereas AF and JL acknowledge 
financial support by the FP7 GTSnext project. 



20 



Appendix 



To illustrate the numerical procedure we have followed to obtain the methods 
with s > 7 stages in this work we now describe in detail the construction of 
method ABAH1064 of generalized order (10,6,4), suitable to be used in Heliocen- 
tric coordinates. 

We consider nine-stage methods with the following sequence of coefficients: 

a\ b\ a 2 b 2 03 63 04 64 05 65 a 5 64 04 63 03 b 2 a 2 h a\ . (32) 

Symmetric splitting methods of generalized order (10,6,4) for Heliocentric co- 
ordinates must satisfy ten order conditions, that is, consistency 

ai + a 2 + a 3 + a 4 + a 5 = -, 2(6i + b 2 + 63 + 64) + b 5 = 1, 



the special constraint (31 ) for Heliocentric coordinates, and the order conditions 
related to the Lyndon multi-indices (3), (5), (7), (9), (1,2), (1,4), (2,3) inTable[l] 
(recall that the q are given by ( 12 )) . We thus have ten polynomial equations and 
ten unknowns, which we collect in a vector x = (ai, . . . , 05, 61, . . . , 65) G M 10 . 
The system of algebraic equations one aims to solve can then be written in 
the compact form f(x) = 0. Recall that any solution of such system must 
have at least one negative cij and one negative bj. We are interested in finding 
solutions with a small Euclidean norm = ||(ai, ■ ■ • , 05, 61, ■ • • , 65) || (so that 
the negative dj and bj have small absolute values). 

In order to do that, we first split the system f(x) = into 

/i(x) = 0, / 2 (x) = 0, (33) 

where f 2 {x) = corresponds to the conditions for the Lyndon multi-indices 
(1,2), (1,4), and (2,3), whereas fi(x) = collects the remaining seven equations. 
We then proceed as follows: 

• We determine the point x° = (a°, . . . , a®, b®, . . . , 65) G M 10 as the (unique) 
solution of the following constrained minimization problem: 



mm 



a 3 =a 4 =0J 1 (a 1 ,...,a 5 ,b 1 ,...,b 5 )=0 . 

It is not difficult to check that the sequence of coefficients 

a? 6? a° 2 (6° + b° + 6°) 4 bl 4 (b° 2 + b° 3 + fig) 4 6? 4 

corresponds precisely to the 5-stage symmetric ABA method of general- 
ized order (10,2) with positive coefficients considered in [TU1 Q2J [T3] . Thus, 
it is our starting point in the search of an efficient (10,6,4) method. 

We then choose an arbitrary orthogonal matrix M £ M 3x10 , and for a 
randomly chosen complex number 7 G C, consider the following one- 
parameter family of systems of polynomial equations: 

/i(a:) = 0, tf 2 (x) + (l-t) 1 M-(x-x°) = 0. (34) 
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For a generic M and 7, there exists a unique continuous curve x = pit) G 
C 10 of solutions of this family of polynomial systems such that p(0) = x°. 
Here t £ [0, 1) denotes the continuation parameter. If x = p(l) € M 10 , 



then x is a real solution of the original system (33). We have applied a 
numerical continuation algorithm to compute such a solution for several 
values of 7 € C, and found two real solutions. The solution x with 
smaller norm ||x|| gives the method ABAH1064 of generalized order (10,6,4) 
displayed in Table [4j Notice the small absolute value of both 03 and 04 in 
the resulting scheme (recall that this solution has been obtained starting 
with x° such that a® = a® = 0). 

We have followed a similar procedure to compute solutions starting with x° 
such that a- 1 = a° = for indices (i, j) 7^ (3,4). Such procedure leads to non- 
real solutions x for some of the choices of and for other choices gives real 
solutions with larger norm ||x|| and larger error terms than method ABAH1064. 
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